module hmdamod
	use myfuncs
	implicit none
	

	integer, parameter	:: nrecs=7
	integer				:: m=16332987,nmsa,ic
	integer, allocatable	:: idata(:,:)
	
	private	:: nrecs,m,nmsa,idata,ic

	contains
	
	subroutine hmdastoredata(filename)
		character(len=*):: filename
		integer	:: ioerr,myunit
        character(len=170) :: line
		
		integer, parameter	:: reclist(2,nrecs)=reshape([82,82,23,23,25,25,43,44,60,60,27,31,86,89],[2,nrecs])	! gender, purpose, owner-occupancy, ST, race
		integer				:: i,j

		open (file=filename,status='old',iostat=ioerr,newunit=myunit)
		if(ioerr>0) then
			print *,"trouble opening file:", filename
			stop
		endif
		m=0
        do
            read (myunit,*,iostat=ioerr)
            if (ioerr/=0) exit     
            m=m+1
		enddo
		print *,"file has ",m," lines"
        rewind(myunit)
		allocate(idata(nrecs,m))
		do i=1,m
            read(myunit,"(a)",iostat=ioerr) line
			do j=1,nrecs
				read(line(reclist(1,j):reclist(2,j)),*,iostat=ioerr) idata(j,i)
				if(ioerr>0) idata(j,i)=huge(1)
			enddo
		enddo	
		close(myunit)
		open(file="c:/dropbox/mystuff/heavymean/hmdamaster.fin",form="unformatted",access="sequential",newunit=myunit,iostat=ioerr)
		if(ioerr>0) then
			print *,"trouble opening file to save"
			stop
		endif
		write(myunit) idata
		close(myunit)
		print *,"done"
	end subroutine

	subroutine hmdaloaddata
		integer	:: ioerr,myunit
		
		open(file="c:/dropbox/mystuff/heavymean/hmdamaster.fin",form="unformatted",access="sequential",newunit=myunit,iostat=ioerr)
		if(ioerr>0) then
			print *,"trouble opening masterfile"
			stop
		endif
		if(allocated(idata)) deallocate(idata)
		allocate(idata(nrecs,m))

		read(myunit) idata
		print *,shape(idata)
		print *,"done loading hmda data"
		print *,"top coded loan amounts",count(idata(6,:)==99999.and.idata(3,:)==1)
		print *,"top coded income",count(idata(7,:)==9999.and.idata(3,:)==1)
		ic=0
		call transcodes
	end subroutine
	
	subroutine transcodes
		integer, allocatable	:: mslist(:)
		integer	:: i,j
		allocate(mslist(0))
		do i=1,m
			if(idata(5,i)<huge(1)) then
				idata(5,i)=merge(1,merge(3,2,idata(5,i)>5),idata(5,i)==5)
			endif
			if(idata(4,i)==huge(1)) cycle
			j=1
			do 
				if(j>size(mslist)) then
					mslist=[mslist,idata(4,i)]
					idata(4,i)=j
					exit
				endif
				if(mslist(j)==idata(4,i)) then
					idata(4,i)=j
					exit
				endif
				j=j+1
			enddo
		enddo
		nmsa=size(mslist)
		print *,"number of regions",nmsa
	end subroutine

	subroutine sethmdadata(vec,randomflag)
		real, allocatable	:: vec(:)
		logical				:: randomflag
		integer, parameter	:: tind=7
		real	:: x(4)
		integer	:: inds(4),i,j,k

		do
			call rnun(x)
			if(randomflag) then
				inds=[2*x(1),3*x(2),2*x(3),nmsa*x(4)]+1
			else
				inds=getindsfromi(ic)			
				ic=ic+1
			endif
			j=0
			do i=1,m
				if(all(idata(1:4,i)==inds(1:4)).and. idata(tind,i)<huge(1)) j=j+1
			enddo
			if(j>=5000) exit
		enddo
		if(allocated(vec)) deallocate(vec)
		allocate(vec(j))
		j=1
		do i=1,m
			if(all(idata(1:4,i)==inds(1:4)) .and. idata(tind,i)<huge(1)) then
				vec(j)=idata(tind,i)
				j=j+1
			endif
		enddo
		print *,"popmax",maxval(vec)
		print *,inds(:)
	end subroutine
	
	subroutine countsubpops
		integer, parameter	:: tind=7
		integer	:: inds(4),c,i,j,i0

		c=0
		do i0=0,2*3*2*nmsa-1
			inds=getindsfromi(i0)
			j=0
			do i=1,m
				if(all(idata(1:4,i)==inds(1:4)).and. idata(tind,i)<huge(1)) j=j+1
				if(j>=5000) exit	
			enddo
			if(j>=5000) c=c+1
		enddo
		print *,"number of subpopulations",c
	end subroutine
	
	function getindsfromi(i0) result(inds)
		integer	:: i0,inds(4),i
		i=i0
		inds(1)=mod(i,2)
		i=i/2
		inds(2)=mod(i,3)
		i=i/3
		inds(3)=mod(i,2)
		i=i/2
		inds(4)=i
		inds=inds+1
	end function

end module